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Abstract We present two Bayesian procedures to infer the interactions and external 
currents in an assembly of stochastic integrate-and-fire neurons from the recording of 
their spiking activity. The first procedure is based on the exact calculation of the most 
likely time courses of the neuron membrane potentials conditioned by the recorded 
spikes, and is exact for a vanishing noise variance and for an instantaneous synap- 
tic integration. The second procedure takes into account the presence of fluctuations 
around the most likely time courses of the potentials, and can deal with moderate 
noise levels. The running time of both procedures is proportional to the number S of 
spikes multiplied by the squared number N of neurons. The algorithms are validated 
on synthetic data generated by networks with known couplings and currents. We also 
reanalyze previously published recordings of the activity of the salamander retina (in- 
cluding from 32 to 40 neurons, and from 65,000 to 170,000 spikes). We study the 
dependence of the inferred interactions on the membrane leaking time; the differences 
and similarities with the classical cross-correlation analysis are discussed. 



1 Introduction 

Over the past decades, multi-electrode recordings (Taketani and Baudry, 2006) have 
unveiled the nature of the activity of populations of neural cells in various systems, 
such as the vertebrate retina (Schnitzer and Meister, 2003), cortical cultures (Tang et 
al, 2008), or the prefrontal cortex (Peyrache et al, 2009). The observation of substantial 
correlations in the firing activities of neurons has raised fundamental issues on their 
functional role (Romo, Hernandez, Zainos and Salinas, 2003; Averbeck, Latham and 
Pouget, 2006). From a structural point of view, a challenging problem is to infer the 
network and the strengths of the functional interactions between the neural cells from 
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the spiking activity (Fig. [T}V). Powerful inference procedures are needed, capable to 
handle massive data sets, with millions of spikes emitted by tens or hundreds of neurons. 

A classical approach to infer functional neural connectivity is through the study of 
pairwise cross-correlations (Perkel, Gerstein and Moore, 1967; Aersten and Gerstein, 
1985). The approach was applied in a variety of neural systems, including the audi- 
tory midbrain of the grassfrog (Epping and Eggermont, 1987), the salamander retina 
(Brivanlou, Warland and Meister, 1998), the primate and rat prefrontal cortex (Con- 
stantidinidis, Franowicz and Goldman-Raking, 2001; Fujisawa, Amarasingham, Har- 
rison and Buzsaki, 2008). Other approaches, capable of taking into account network- 
mediated effects, were proposed based on concepts issued from statistics and graph the- 
ory (Seth and Edelman, 2007; Dahlhaus, Eichler and Sandkiihler, 1997; Sameshima and 
Baccala, 1999; Jung, Nam and Lee, 2010), information theory (Bettencourt, Stephens, 
Ham and Gross, 2007), or statistical physics (Schneidman, Berry, Segev and Bialek, 
2006; Shlens et al, 2006). 

An alternative approach is to assume a particular dynamical model for the spike 
generation. The generalized linear model, which represents the generation of spikes as 
a Poisson process with a time-dependent rate is a popular framework (Brown, Nguyen, 
Frank, Wilson and Solo, 2001; Truccolo et al, 2005; Pillow et al, 2008). The Integrate- 
and-Fire (IF) model, where spikes are emitted according to the dynamics of the mem- 
brane potential is another natural candidate (Gerstner and Kistler, 2002; Jolivet, Lewis 
and Gerstner, 2004). The problem of estimating the model parameters (external cur- 
rent, variance of the noise, capacitance and conductance of the membrane) of a single 
stochastic IF neuron from the observation of a spike train has received a lot of at- 
tention (Paninski, Pillow and Simoncelli, 2004; Pillow et al., 2005; Mullowney and 
Iyengar, 2008; Lansky and Ditlevsen, 2008). Few studies have focused on the inference 
of interactions in an assembly of IF neurons (Makarov, Panetsos and de Feo, 2005). 
Recently, we proposed a Bayesian algorithm to infer the interactions in a network of 
stochastic perfect integrators when the synaptic integration is instantaneous and the 
noise is vanishingly small (Cocco, Leibler and Monasson, 2009). 

In the present work we introduce a Bayesian algorithm to infer the couplings and 
the external currents in an assembly of leaky IF neurons, and in presence of moderate 
input noise (Fig. [T}\). The computational time grows as the product of the number of 
recorded spikes, and the square of the number of neurons. We validate the algorithm 
on synthetic data, and apply it to real recordings of the ganglion cell activity in the 
salamander retina, presented with natural visual stimuli, and in the absence of stimulus 
(spontaneous activity). 

2 Materials and Methods 

2.1 Definition of the Leaky Integrate-and-Fire model 

In the Leaky Integrate-and-Fire (LIF) model, the membrane potential Vi(t) of neuron 
i at time t obeys the first-order differential equation, 

c^(t) = - 9 Vi(t) + i° yn (t) + h + m(t) (i) 

where C and g are, respectively, the capacitance and conductance of the membrane. 
The ratio r = C/g is the membrane leaking time. I^ yn (t) is the synaptic current coming 
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Fig. 1 A. Extra-cellular recordings give access, through spike-sorting, to the times t, ^ of the 
spikes emitted by a population of neurons. We want to infer the values of the interactions 
Jij and external inputs /; of the network most likely to have generated the recorded activity. 
B. Example of firing activity of N = 2 neurons. The top panel shows three spikes emitted by 
neuron 1. Panels a, b, c show possible activities of neuron 2, with equal average firing rates 
but with different timings. 



from the other neurons and entering the neuron i at time t: 

irw= £ - y « £*(*-**,*) ( 2 ) 

where is the strength of the connection from neuron j onto neuron i (Figure UJA.); 
tj i- is the time at which neuron j fires its k th spike. We assume that synaptic inputs are 
instantaneously integrated, i.e. that the synaptic integration time is much smaller than 
all the other time scales, including r. Our method for inferring the interactions relies 
on this assumption, and should be modified in the presence of synaptic integration 
kernels with temporal filtering, /j is a constant external current flowing into neuron 
i (Fig. [TJ\), and r)i(t) is a fluctuating current, modeled as a Gaussian noise process: 
(Viit)) = 0, {ni(t) Vj(t')) = °" 2 $ij — 0- The noise standard deviation, a, has here 
the dimension of a current times the square root of a time. An alternative definition 
would consist in rescaling a with a time dependent factor, e.g. \fr; our definition allows 
us to reach the perfect integrator limit (r — > oo) while keeping a fixed. 

The neuron i remains silent as long as Vi remains below a threshold potential Vj/j. 
If the threshold is crossed at some time to, i.e. Vj(trj) = Vth, then a spike is emitted, 
and the potential is reset to its resting value: V(£q") = 0. The dynamics then resumes 
following ([TJ. 



2.2 Likelihood of the spiking times for given interactions and currents 

Let J = {Jij} and X = {/j} denote the sets of, respectively, the interactions and 
currents. Let f. G [0; T] be the time at which neuron i emits its k th spike; T is the 
duration of the recording. How can we infer the interactions and currents from the 
observation of the spiking activity? Consider the raster plots in Fig. [TJ3. In pattern a, 
the timings of the spikes of neuron 1 do not seem to be correlated to the activity of 
neuron 2. Hence, we may guess that there is no interaction from neuron 2 to neuron 
1 (Jl2 = 0). In pattern b, a spike of neuron 1 is likely to follow a spike of neuron 2, 
which suggests that the interaction J12 is positive. Conversely, in pattern c, it seems 
that the firing of neuron 2 hinders the firing of neuron 1, which indicates that J12 has 
a negative value. 
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This crude reasoning can be made mathematically rigorous in the framework of 
statistical inference (Cover and Thomas, 2006). Let us define the likelihood P(T\J ,X) 
of a set of spiking times, T = {t,; given J and X. According to Bayes rule the most 
likely couplings and currents, J and X, given the set of spiking times T can be inferred 
through the maximization of P(T\J ,X) Q. Due to the statistical independence of the 
noises 77^ from neuron to neuron, the likelihood P of T given J ,X can be written as 
the product of First-Passage Time (FPT) probabilities, 

P(T\J,X) = WvFPTik.k+iKk, {t j4 }, {Jij}, h) . (3) 

i,k 

Here Pfpt denotes the probability that Vi crosses Vth for the first time at time k+l> 
starting from at time % and conditioned to the inputs from the other neurons at 
times tj g, with t^k < tj i < tj fc+i. As the synaptic integration is instantaneous, an 
incoming spike from neuron j results in a (positive or negative) jump of the potential Vi 
by Jij/C; Pfpt can therefore be interpreted as the FPT density probability for a one- 
dimensional Ornstein-Uhlenbeck process with a time-dependent force. It is important 
to stress that the presence of the products over the spike intervals in Q does not entail 
that the spiking times are independent. 

Consider now the potential Vi(t) during the inter-spike interval (ISI) [ijfejtj fc+l]- 
The boundary conditions are Vi(t^ k ) = (reset of the potential right after a spike), 
and Vi(t~ k+1 ) — Vth (condition for firing). At intermediate times, the potential can 
take any value smaller than Vth- The logarithm of the probability of a dynamical path 
(time course) of the potential over the k th ISI of neuron i is, after multiplication by 
the variance a 2 of the noise, 



C[Vi(t);k,T,J,X\ = -\ f Zk+1 dt m (tf (4) 

J ti,k 



C^(t)+g Vl (t)-ir(t)-h 2 



according to the Gaussian nature of the noise %(£) and to the dynamical equation of 
the LIF ©. 



2.3 Dynamical equations for the optimal potential and noise 

While no exact expression is known for Pfpt, h can be analytically approximated 
by the contribution of the most probable dynamical path for the potential, V*(t) 
(Paninski, 2006). This approximation becomes exact when the standard deviation a 
of the noise is small. The idea is to replace the distribution of paths for the potential 
Vi(t) with a single, most likely path V*(t), which we call optimal. We now explain how 
to derive V*(t) through the condition that the log-probability C ((4]) is maximal. 

Let us assume first that V*(t) < V t h- Then, the derivative of C in Q with respect 
to V*(t) must vanish, which gives 

■2 T/ * HT syn 
= - C -tt^ + 9 2 V * (*) + - 9 in*) ~9li = 0. (5) 



8C 
6Vi(t) 



We consider here that the a priori measure over the couplings and currents is flat. 
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We now turn this second order differential equation for the optimal potential into a 
first order differential equation at the price of introducing a new function, rj*(t), and 
a new first order differential equation for this function. It is straightforward to check 
that the solution of 

HV* 

C M (t) = ~ 9 V ' {t) + T i W + h + V:{t) (6) 
is a solution of the optimization equation (J5]l if r\l (t) fulfills 

M W = ^ ?W = ^), (7) 

where r is the membrane leaking time. The similarity between eqns {TJ and © allows 
us to interpret r/*(t) as a current noise. However, this noise is no longer stochastic, 
but rather it follows the deterministic path solution of iff}. We will, therefore, in the 
following refer to r]* (t) as the optimal noise, rj* (t) corresponds to the most likely value 
the noise takes given the set of spiking times. Solving ([7| shows that the optimal noise 
is an exponential function of the time: 

rfiit) = V exp(+t/r) {V*{t) < V th ) , (8) 

where 77 is a constant, which we call noise coefficient. 

It may happen that the optimal potential only reaches the threshold without ac- 
tually crossing it at intermediate times. When this is the case, the optimal potential 
equals V* (t) — V t ^ and its derivative with respect to the time vanishes. The expression 
for the optimal noise can be then read from (|SJ, and is given by 

v* (t) = gv th - it vn {t) - h (v*(t) = v th ) • (9) 

Equation © ensures that the potential does not cross the threshold value at a time 
t < ti,k+l- 

Despite their apparent simplicity, eqns H6I8I9[) ar e not easy to solve, due mainly 
to the interplay between the two regimes, V* < V t h an d V* — V t h, mentioned above. 
The determination of V*(t) was achieved numerically by Paninski for a single neuron 
(Paninski, 2006). We now sketch the procedure to determine V*(t) rapidly, even for 
tens of neurons. The procedure relies on the search for contacts, that is, times at which 
the optimal potential touches the threshold. There are two types of contacts: contacts 
coinciding with a synaptic input (the potential touches the threshold at time tjj.), 
and contacts arising in between two inputs. In the absence of leakage, only the former 
type of contacts matter, and a search procedure to locate those isolated-time contacts 
was proposed by Cocco, Leibler and Monasson (2009). In the presence of leakage, both 
types of contacts have to be taken into account. The search procedure is more complex, 
and is explained below. 

2.4 Fixed Threshold procedure: optimal paths for the potential and the noise 

We assume in this Section that the couplings and currents are known. Consider neuron 
i at time t £ fcjtj fc+i], where k is the index of the ISI. The initial and final con- 
ditions for the optimal potential are: V*(t~^ k ) = and V* (t~ fe+1 ) = Vth- in between, 
V* (t) obeys the LIF evolution equation ([6]) with an optimal 'noise' rj*(t). r)*(t) can 
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Fig. 2 A & B. Sketches of the optimal potentials V* (top, black curves) and noises 77* 
(bottom, black curves) for one neuron receiving several inputs from two other neurons (red 
and green impulses, middle, with J re( j = — J green = .2CVth). The membrane conductance 
is g = .81 /Vth- The jump in the optimal noise consecutive to an active contact is always 
positive. A. Illustration of Passive (P) and Active (A) contacts. B. Comparison with numerical 
simulations, averaged over ~ 5,000 samples, for <r = .07 (red), .18 (purple), .36 (blue); noise 
curves are averaged over the time-window At = .015 T. C. Dashed lines represent possible 
paths for the potential when the noise standard deviation, <r, docs not vanish. The amplitude 
of the fluctuations of the potential around V at the mid-point of the ISI is symbolized by 
the double arrow line. D. Probability p B (St\V) that an Ornstein-Uhlenbeck process starting 
from V does not cross the threshold Vth f° r a time <5f = 3.5 r. Parameters are gV t h/I = 1.2, 
<t = .15. The tangent line to p s in V = Vth crosses the p a = ^ l me m V-th ■ System of two 
IF neurons, with gV th /h = l-5,gV th /I 2 = 2., J 12 /(CV th ) = .1, J 2 i = 0,5" = .25. The dashed 
and full black curves represent the optimal potentials for neuron 1 calculated by, respectively, 
the Fixed and Moving Threshold procedures; for the latter, Vju is shown in red. One random 
realization of the membrane potential (averaged over a 10 msec time-window) is shown for 
comparison (blue curve). 



be interpreted as a non-stochastic, external, time-dependent current to be fed into the 
neuron in order to drive its potential from to V t h, given the synaptic couplings. The 
expressions for the optimal 'noise' are given by (JH]) when V*(t) < Vth, an d © when 
the optimal potential (t) is equal to the threshold value. 

When V*(t) reaches the threshold at a time coinciding with an incoming spike, 
the coefficient 77 in (|8l) may abruptly change through an active contact; the notion of 
active contact is illustrated in the simple case of a neuron receiving a single spike in 
Appendix lA.il The potential V*(t) may also touch the threshold without crossing it, 
and the noise may remain constant over some time interval; we call such an event 
passive contact. That the potential can brush, or remain at the threshold level without 
producing a spike is made possible by the a — > limit. We will discuss later on the 
validity of this calculation, and how to modify it when the noise standard deviation, 
a, does not vanish. Both types of contacts are shown in Fig. 

Let us explain how the positions of active and passive contacts can be determined. 
Let t\ < t2 < ■ ■ ■ < t^j be the emission times of the spikes arriving from the neurons 
interacting with i during the time interval [to = £j fc; ^M+l = ij fe+l]> an d J\ , J2, ■ ■ ■ , Jm 
the corresponding synaptic strength^. Let Vq — be the initial value of the potential, 
and mo = 1 be the index of the first input spike. If the time is small enough the optimal 
potential is surely below the threshold value. According to ([§} the optimal noise is an 

2 Due to the limited temporal resolution of the measurement two inputs of amplitudes J 
and J' can apparently arrive at the same time; if so, we consider, based on model and 
that a single input of amplitude J + J' enters the neuron. 
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exponential with noise coefficient rj, and the optimal potential is obtained by solving 
(O with the result, 

M 

V(v, t) = V e-C-^ + £ ^ e -(*-*™)/^( 4 - t rn ) 

m— mo 

+ £i ( i_ e -(*-M/T ) + i? sinh (^M (10) 

where is the Heaviside function. It is tempting to look for the value of n such that 
a spike is emitted at time tM+i, defined by the implicit equation Vi(r), tjtf+l) = Vth- 
However, the corresponding potential might not be below threshold at all intermediate 
times to < t < tM+l- Instead, we look for the smallest noise capable of driving the 
potential from its initial value V$ into contact with the threshold: 

?7* = min{?7: max Vj(jj, t) = V t h} ■ (H) 

to<t<tM + l 

As the potential (|10[) is a monotonically increasing function of rj, a value of the noise 
smaller than rf would not be able to bring the potential to the threshold and to 
trigger a spike at any time, while a value larger than rf would violate the condition 
that the potential cannot cross the threshold on the time interval to < t < t^+i, see 
Appendices IA.2I and IA.3I 

We denote by t c the time at which the threshold is reached: Vi(rf ,t c ) = Vth- 
The solution to the minimization problem (|ll|l can be found following the procedure 
described below. Briefly speaking, the procedure identifies candidates for the contact 
points, selects the best one, and is iterated until the ISI is completed. 

— Active candidates: we first consider the possibility that the contact time t c coincides 
with a synaptic input. We therefore calculate for each m = mo, . . . , M + 1, the 
root ?7 m of the implicit equation V(r) m ,trn) = Vth- The smallest of those M noise 
coefficients is called rf. 

— Passive candidates: we then consider the case where the contact time t c may not be 
simultaneous to any input, but rather fall between two successive spikes. For each 
< m < M, we look for a noise coefficient r/ p and a contact time t c € \t m ',tm+l\ 
fulfilling the set of coupled equations Vi(rj p ,t c ) = Vth,Vi(rip,t c ) = 0, expressing 
that the potential reaches and does not cross the threshold. These two equations 
can be solved analytically, see expressions (|51|l and (|52p in Appendix [B] We call r]p 
the smallest noise coefficient corresponding to those possible passive contacts. 

— Selection of the best candidate: 

• if rja < r]p, the contact is active and takes place at time t c = t m * for a 
certain m* comprised between tuq and M + l (Fig. UK). The optimal potential and 
noise in the time interval [to, t m *\ are given by, respectively, eqns ()10|) and |(5J with 

v = n ■ 

• If r]p < rf, the contact is passive, and takes place in the time interval 
\bm c -V,tmc\ for a certain m c comprised between mo and M + 1 . The potential 
will remain equal the threshold, and the noise will remain constant according to 

over a finite time interval [t c ;t c + A c ], after which both V*(t) and rf{t) resume 
their course (Fig.[2]4.). A c is the smallest delay allowing the potential to be in active 
contact with the threshold at a later time t m *, with m* > m c . The correctness 
of this statement is ensured by the fact that there can be at most one passive 
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contact between two active contacts (Appendix [EJ; hence, a passive contact is 
necessarily followed by an active contact (or by the spike at the end of the ISI). 
For every integer m comprised between m c and M + 1, we calculate analytically 
the delay A c (m) such that the potential reaches the threshold in t m , see eqn (|53|l 
in Appendix [B] the smallest among those delays and the corresponding value of m 
are, respectively, A c and m*. 
— Iteration: we are left with the calculation of r\* (t) and V* (t) on the remaining part 
of the inter-spike interval, [t m * \ tM+l]- To do so, we iterate the previous steps. 
We first update to «- t m «, m <s— m* + 1, Vo *- V th + 6(—J m »)^- in (JTOj) , and 
look for the lowest noise producing a new contact over the interval [to, tjvf+l] using 
again. The procedure is repeated until the whole inter-spike time interval is 
exhausted. 

As a result a sequence of values for rf is built, each value corresponding to the noise 
coefficient (|lip between two successive active contact points. 



2.5 How small should the variance of the noise be? 

The L1F dynamical equation fTJ) involves quantities, such as the membrane potential, 
the membrane conductance, the input current, which have different physical units. 
A straightforward algebra shows that {TJ is equivalent to the following differential 
equation, 

% = -v + E J v E s $- hk) + + nm , (12) 

which involves only dimensionless variables (denoted with overbars): t = t -jj, V = 

Y^i Jij — cVth ' ^ = g Vth ' ^ ne n ™ se ^ as zero mean, and covariance {fji(t)fjj(P)) = 
a 2 5ij8(i — i'), where 

(13) 



VthVgC 

Intuitively, we expect that the potential Vi(t) will not depart much from the optimal 
path V*(t), and, hence, that our inference algorithm will be accurate if the dimen- 
sionless standard deviation of the noise, a, is small. We illustrate this claim on the 
simple case of a neuron receiving a few inputs from two other neurons during two 
inter-spike intervals (ISI) of length 2r, see Fig. fjj3. The times of the input spikes were 
randomly chosen, once for all, before the simulations started. Then, we numerically in- 
tegrated the LIF equation for the potential |T]) for 10 6 random realizations of the noise 
rj(t). The realizations such that the neuron spiked twice, with ISIs falling in the range 
[1.99, 2.01] x r were considered as successful. The number of successful realizations was 
comprised between 10 3 and 10 4 , depending on the noise level, a. We show in Fig. [2j3 
the paths of the potential and of the noise, averaged over successful realizations, and 
compare them to the optimal potential, V* , and noise, rf . As expected the agreement 
is very good for small a. We now make this observation quantitative. 

Consider the k th inter-spike interval \pi k>ti k+l] °f neuron i. The optimal potential 
V*(t) is the time-course followed by the LIF membrane potential Vi(t) in the a — > 
limit. When the noise variance is not vanishing, the potential Vi(t) can slightly deviate 
from the optimal path (Fig.[5p). Deviations are null at the extremities of the inter-spike 
interval due to the boundary constraints on the potential. A measure of the magnitude 
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of the fluctuations of the potential is thus given by the variance of Vi(t) — V*(t) at the 
middle of the ISI, i.e. t = ^(t i>k + t i;k+1 ) (Fig.|2p). This variance can be calculated 
when the constraint that the fluctuating potential Vi(t) does not cross the threshold 
at times t < tj is relaxed, see Appendix iDl We obtain 

w-v?f}_ & 2 tanh fw-^y (14) 



where r is the membrane leaking time. As expected, if a is small, so are the fluctuations 
of the potential around the optimal path. 

However, the reverse statement is false. Consider, for instance, the case of a perfect 
integrator, for which the dimensionless a (|13|) is virtually infinite. Sending g — > in 
(flUl , we obtain 

vi 2{cv th y {9 ^ 0) - (15) 

Hence, the relative fluctuations of the potential are small if the typical amplitude of the 
electrical charge entering the neuron during the ISI due to the noise, cr^/i^fc+i — tj fc, 
is small compared to the total charge C Vth necessary to reach the threshold from the 
rest state. It is interesting to note that this statement applies to the LIF, too. Whatever 
the level of the noise, a, the relative fluctuations of the potential (|14p can be made 
small if the duration of the ISI is short enough compared to the membrane leaking 
time, r. 



2.6 Beyond the weak-noise limit: the Moving Threshold procedure 

For large values of a, a discrepancy between the optimal potential and the potential 
obtained in simulations appears (Fig. 03). A general observation is that the optimal 
potential calculated by the Fixed Threshold procedure can get very close to Vth, while 
the true potential stays further away from the threshold to avoid premature firing. To 
further illustrate this effect, consider a system of two IF neurons, 1 and 2, both fed 
with an external current. In addition, neuron 1 receives positive inputs from neuron 
2 (J12 > 0), and neuron 2 is independent from the activity of neuron 1 (J21 = 0). In 
presence of a strong noise, the optimal potential calculated from the Fixed Threshold 
procedure quickly reaches a stationary value close to V t h, while the random potential 
obtained from simulations fluctuates around a much lower level (Fig. 03). The presence 
of a strong noise biases the membrane potential to lower values to prevent early spiking. 
A heuristic approach to reproduce this bias consists in decreasing the threshold from 
V^ to a time- and context-dependent value, V t ^f . We now explain how this moving 
threshold, VM , is determined. 

Consider first a neuron with no synaptic input, fed with an external current I, 
during the inter-spike interval [t{ k ; k+l]- We call p s (5t\V) the probability that the 
potential, taking value V at time fc+1 — remains below the threshold at any larger 
time t, with t i — St < t < fe+i ■ This probability depends on the current /, and 
can be expressed for an arbitrary level of noise, a, as a series of parabolic cylinder 
functions (Alili, Patie and Pedersen, 2005). Figure [2p show p s as a function of V for 
some characteristic values of the parameters. The probability of survival, p s , sharply 
decreases to zero when V gets close to the threshold, V — V t h- We model this outcome 
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by the following approximation, which involves a new, effective threshold V th : we 
consider that the processes starting from a value of the potential V > Vjjf will not 
survive for a time delay St. In other words, the true threshold, V t h, is changed into a 
'moving' threshold, which is a function of the current 7, the time St, and the parameters 
g, C, a. A simple way to define VM is to look at the intersection of the tangent line 
to p s in V = Vth with, say, the p s = \ linfl the resulting expression for V^jf is given 
in Appendix [Ej Figure 03 shows the output of the Moving Threshold procedure on 
the simple 2-neuron system described above. The optimal potential, 'pushed' down by 
the moving threshold Vjjf is much lower than in the Fixed Threshold approach and in 
much better agreement with the random realization of the membrane potential. More 
details are given in Section r3.1.4l 

To account for the existence of synaptic inputs, we may choose the parameter 
/ entering the calculation of p s and V t ^f to be the value of the effective current 
If = h + Jij fj > rather than the external current Jj itself. Here, fj is the av- 

erage firing rate, defined as the number of spikes fired by neuron j divided by the 
duration T. Contrary to the external current 7j, the effective current If takes into 
account the (average) input current coming from other neurons. This choice was done 
in the numerical experiments reported in the Results section. To further speed up the 
calculations, we derive the value of V^Jf for discrete-values delays St only; in a discrete 
interval, Vjf is kept to a constant value. 

Alternative heuristic approaches to deal with the presence of moderate noise can be 
proposed. In Appendix[E] we introduce a cost-function for the effective current, whose 
effect is also to decrease the optimal potential. These approaches are effective when 
the optimal potential calculated by the Fixed Threshold procedure quickly saturates 
to a level close to V t h- More precisely, we expect the Moving Threshold procedure to 
be efficient if the membrane leaking time is smaller or comparable to the ISI, and the 
leaking current, ~ gV t h, is larger or equal to the external current, I. 



2.7 Maximization of the log-likelihood to infer the interactions and currents 

The Fixed or Moving Threshold procedures allow us to calculate the optimal paths 
for the potential and the noise, given the couplings and currents. Knowledge of those 
paths gives us also access to the logarithm of the likelihood P in the a — » limit, 

L*(T\J,X) = lim a 2 log P(T\ J, I) 

<T->0 

= $>[t?(t);*,r,.7,Z] = -|X; [ dtr,*(tf (16) 

i.k i.k •' ti - k 

Since L* in (1161) involves the sum over different neurons, the maximization over the 
couplings Jij and the current Ii of neuron i can be done independently of the other 
couplings Jj/j and currents JV (i 7^ i). Formally, we are left with Af independent 
inferences of the most likely couplings and current for a single neuron, in presence of 
the spikes emitted by the A^ — 1 other neurons. As a consequence neurons 'decouple' in 

3 This choice is arbitrary; other values, ranging from j to 1 have been tried, do not quali- 
tatively affect the results presented later in this article. 
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the inverse problem: the couplings Jij and the current Jj of neuron i can be inferred 
independently of the other couplings J^j and currents 1^ (i 7^ i). 

C defined in Q is a negative-semidefinite quadratic function of its arguments 
Vi(t), Jij , Ii. It is thus a concave function of the couplings and the currents. This prop- 
erty holds for L* (|16|l (Boyd and Vandenberghe, 2004). In order to infer the most likely 
current Ii and couplings Jy, we start from an arbitrary initial value e.g. Ii — Jij = 0. 
The full path of the optimal noise, rj*(t), over all the inter-spike intervals k of neuron 
i, is calculated following the above procedure. We then update the couplings and the 
current using the Newton- Raphson method to maximize logP, i.e. to minimize the 
integral of the squared optimal noise, see H16[) . Convergence follows from the concavity 
property stated above. The procedure requires the expressions for the gradient and the 
Hessian matrix of logP with respect to the couplings Jij and the current Ii, which 
can be calculated exactly from (|16p and JTUJ). Note that logP is piecewise continuously 
twice-differentiable; while the gradient is continuous for all Jij and Ii, the Hessian 
matrix is bounded and negative, and may discontinuously jump due to a change of the 
contact points. Knowledge of the Hessian matrix is also important to determine how 
reliable are the values of the inferred parameters. 



2.8 Accuracy on the inferred parameters 

When the variance of the noise, a 2 , vanishes the inferred parameters cannot deviate 
from their most likely values. However, for small but non zero a, deviations are possi- 
The probability for such deviations can be estimated from the expansion of L* 
around its maximum. We introduce for each neuron i, the iV-dimensional vector w, 
whose components are: = Zj r, and Vj = Jij for j 7^ i. The multiplication of the 
current by the membrane leaking time ensures that all components can be expressed 
in units of a coupling. Similarly we call the vector obtained when the current and 
couplings take their most likely values, that is, maximize L*. Let us call 

the Hessian matrix of L . The parameters v]j are normally distributed around their 
most likely values, with a covariance matrix given by 

((«f-«f)(«i?-^))=[HW]7j,. (18) 

In particular, the error bars on the inferred parameters are given by the diagonal 
elements of the inverse of H*- 1 -*. Note that, if the value of a is not known, formulas (|17l) 
and (|18[) can still be used to compare the error bars between each other. 

As the entries of scale linearly with the duration T of the recording, or, more 
precisely, the number S of recorded spikes the uncertainty on the inferred parameters 
will decrease as S 1-1 / 2 . A detailed spectral analysis of tj 2 H^ /S in the case of weak 



4 Note that the inferred parameters might be less sensitive than the time course of the 
potential to the noise level a. The reason is that the corrections to the log-likelihood L* , 
to the lowest order in the noise variance a 2 , do not depend on the current and interactions 
(Appendix [D}. 
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couplings, reported in Appendix[Cl shows that the largest eigenvalue, Xmax, is related 
to the fluctuations of the effective current, 

If I, I Y. J » ff > M! " 

when- 

/;• J £ exp(-^±±-^) (20) 

is the average firing rate of neuron j, calculated over the time scale ~ min(r, I SI) 
preceding a spike of neuron i. The smallest eigenvalue, A m j„, corresponds to the fluc- 
tuations of the current 7j alone. In other words, the uncertainty on the inferred value 
for If is much smaller than the one on the current I{. The intermediate eigenmodes 
describe correlated fluctuations of the couplings. Explicit expressions for the largest 
and smallest eigenvalues, Xmax and X m i n , are derived in Appendix [Cl 

When a small change of J and X causes a modification of the set of contact points 
the second derivative of L* may be discontinuous. A simple illustration is provided by 
the the case of a single input, whose log-likelihood L* is reported in Appendix lA.il If 
the maximum is located at, or very close to the boundary dividing two or more sets of 
contacts, the value of the Hessian matrix will depend on the direction along which the 
maximum J, I is approached. This phenomenon is also encountered in the analysis of 
real data, see Section 15.2.21 



3 Results 

3.1 Tests on simulated data 

In this Section, we test our inference procedure on synthetic data generated from 
networks with known interactions and currents. We compare the results obtained from 
our two inference algorithms, the Fixed and Moving Threshold procedures, respectively 
defined in Sections 12.41 and 12.61 

3.1.1 Scaling of the computational time 

We first consider N (ranging from 20 to 160) neurons, with no leakage (g — 0). The 
neurons are uncoupled (J.y = for i ^ j), and fed with identical currents (1$ — I 
for all i). The choice of the noise variance, a 2 , is specified later. The LIF equation is 
solved numerically, using a fourth-order Runge-Kutta integration scheme. We choose 
the elementary time step to be 10 -5 sec, while the average duration of the ISI is 10 3 to 
10 5 longer. For each realization of the noise, the simulation is run until a set of ~ 10 7 
spikes is generated. We then use the first S spikes in this set to infer the currents 
and the couplings (not fixed to zero a priori) with the Fixed Threshold procedure. 
The algorithm stops if the log-likelihood L* increases by less than e = 10~ 12 after 
an iteration of the Newton-Raphson procedure. Alternatively, the algorithm may halt 
when the overall change in the couplings and current becomes smaller than a certain 
a priori bound. 

Figures [3J\&B show how the running time scales with, respectively, the number S 
of spikes, and the number N of neurons. The empirically found scaling, 0(S N 2 ), can 
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Fig. 3 Results of the Fixed Threshold algorithm on a network of TV uncoupled neurons, and 
in the absence of leakage. The running time, on one core of a 2.8 GHz Intel Core 2 Quad 
desktop computer, is shown as a function of the number of spikes, S (A), and of the number 
of neurons, TV (B). C. Inference errors e s on the currents li, If, and on the couplings Jij vs. 
S/N, for TV = 40 neurons and three values of the noise ratio (|22j: r = .4 (A), .04 (■), .004 
(•). Data are shown for one randomly drawn sample; samplc-to-samplc fluctuations are of the 
order of the symbol size. Full lines show square root, linear and quadratic increases (in log-log 
scale); dotted lines serve as guides to the eye. 



be understood as follows. Consider one neuron, say, i. The number of spikes of neuron 
i is, on average, equal to S/N ~ / T, where T is the duration of the recording and / 
is the average firing rate. The number of contact points, N co , is found to scale as the 
number of spikes, S/N . The calculation of the contribution to the Hessian coming 
from the interval between two successive contact points of V* takes 0(N 2 ) time. The 
total calculation of thus requires N co TV 2 ~ S N operation^]. The loop over the 
neuron index, i, gives an extra (multiplicative) factor N. 

The running time of the Moving Threshold algorithm grows as S N 2 , too. However 
the proportionality constant is generally higher than for the Fixed Threshold procedure, 
due to the extra computational burden to calculate Vjjf . For fixed N and S, the running 
times of both procedures increase with the number of contacts, e.g. when the membrane 
conductance g increases. This effect is described in Section [3.21 



3.1.2 Dependence of the inference error on the number of spikes 



We define the inference errors as the root mean square of the difference between the 
inferred parameters, J" 1 ^ , i,- and the true values, J,;, = 0, L = I: 



N(N — 1) 



jinf 



cv th 



1 / I ini 
N^ \ ~T~ 



(21) 



5 Note that the ratio of the time to calculate HW over the time required for the inversion 
of the Hessian matrix is equal to N co N 2 /N 3 ~ S/N 2 , and is generally much larger than one. 
The reason is that the number of parameters to be inferred, TV, has to be smaller than the 
number of constraints over the optimal potential, TV co . For the real data analyzed in Section 
13.21 we have S/N 2 ~ 64 and 108 for, respectively, Dark and Natural Movie data sets. 
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together with a similar definition for the effective current, e(If), with J? replaced 
with the inferred value for If. The inference errors depend on the dimensionless noise 
raticH 

r = a . (22) 

Figure [2P shows the inference errors found for different noise ratios r, and their de- 
pendence on the number S of spikes, in the absence of membrane leakage. For small 
data sets, the inference error is mainly due to the imperfect sampling. As the number 
S of spikes increases, e s decreases as S -1 / 2 , as expected from Section 12.81 When S 
is very large, the errors saturate to a residual value, too- The presence of the residual 
error too results from the dominant-path approximation done in our calculation of the 
likelihood P. The value of too decreases with r as expected. 

The cross-over between the sampling-dominated and residual error regimes takes 
place for a certain value of the number of spikes, S c .o.- Both S c .o. and too depend on 
the observable, i. e. I, I e , J, and on the noise ratio r. With the values of S reached in the 
simulations, the onset of the cross-over is clearly visible for I e , can be guessed for 7, and 
is not observable for J. The existence of a cross-over, and an estimate of S c .o. can be 
derived from the discussion of Section[2]8] When S is large, the a posteriori distribution 
of the inferred parameter, v = /, I e , or J, becomes Gaussian, with a variance 

((Avf) ~ ^1 , (23) 

where A is the eigenvalue of the Hessian matrix of L* attached to the fluctuations 
of the parameter v. The inference error sums up contributions coming from both the 
sampling fluctuations and the residual error. The cross-over takes place when both 
contributions are comparable, ((Av) 2 ) 3 = too, that is, for 

Sc.o. ~ -^r . (24) 

Figure[3p confirms that Sc.o. diminishes with a (or r), and is much smaller for I e than 
for I (as expected from the dependence on the eigenvalue A); moreover, the residual 
error on the couplings is extremely small (or might be even zero). 

As a conclusion, our inference algorithm is very accurate in the absence of mem- 
brane leakage. With 10 3 spikes per neuron only and r = .004, for instance, the errors 
on the currents and on the couplings are, respectively, e s = 3 10 -3 and 4 10 -4 . Even in 
the presence of strong noise (r = .4), and with the same number of spikes per neuron, 
the errors on the effective currents and on the couplings are less than 1%. 



3.1.3 Performance of the Fixed Threshold procedure on networks of coupled neurons 

We now study the ability of the algorithm to infer the interactions between coupled 
neurons. To do so, we consider random connection graphs built in the following way 
(Bollobas, 2001). We start from a complete oriented graph over N neurons, and erase 
each one of the N(N — 1) link with probability 1 — p, independently of each other. 
The removal process is not symmetric: the link i — > j may be removed, while the 

6 When g = 0, changing the value of the current I amounts to changing the time-scale of 
the evolution of the potential in Hence, the errors e s depend on the parameters /, C, o~, Vth 
through the value of r only (as long as / > 0). 
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Fig. 4 Results from the Fixed (empty squares) and Moving (full squares) Threshold algo- 
rithms on a random network of N = 40 coupled neurons; the maximal amplitude of synapses 
is Jo = .2CVth- The error on the couplings, e s (J)/Jo, is plotted as a function of the fraction 
p of connections (A) and conductance over current ratio, gVth/I (B). Each simulated data 
set contains S = 5 10 5 spikes, which is larger than the cross-over size S c . .', the symbol size 
correspond to the fluctuations estimated from ten different data sets for the same network of in- 
teractions. C. Inferred interactions vs. true values of for various values of gVth/I and r, and 
for one random network with a fraction p = .2 of connections. Dashed lines have slope unity. 
Panels C-l, C-2, C-3 show the results of the Fixed Threshold (FT) procedure; the slopes of the 
best linear fits (full lines) arc indicated between parenthesis. Panel C-4 shows the outcome of 
the Moving Threshold (MT) procedure; even if multiplied by 10, the FT couplings of Panel C-3 
are in much worse agreement with the true interactions than the MT couplings. D. Optimal 
potentials V* obtained with the Fixed Threshold procedure for g = .11/Vth (dashed curve) 
and g = 1.21 /Vth (full curve), and for one arbitrarily chosen neuron among the N = 40 neural 
cells; the noise ratio is r = .15. E. Comparison of a random realization of the potential V 
(red) with the optimal potential V (black) obtained with the Moving Threshold VM (green) 
procedure. The network of interactions, the spiking times, and the arbitrarily chosen neuron 
are the same as the ones in D for g = 1.21 /Vth- The time-average of Vjjr is ~ .93 V t h- 



connection j — ^ i is preserved. At the end of the construction process, the average 
number of outgoing (or incoming) neighbors of a neuron is p(N — 1). Each existing 
connection is then assigned a synaptic weight, uniformly at random over the interval 
[—Jo; Jo]- All neurons receive the same external current I. In addition, the membrane 
conductance, g, is now different from zero. The values of p, Jo,I,g, and a are chosen 
so that the network remains below saturation. 

We have also performed simulations where the interaction graph is drawn as above, 
but each neuron i is chosen to be either excitatory or inhibitory with equal probabilities. 
The outgoing interactions from i have all the same sign, and random amplitudes in 
[0; Jo]. The performance of our inference algorithms are qualitatively similar for both 
models. 

Figure |4j\ shows the error on the couplings inferred with the Fixed Threshold 
algorithm, e s (J), as a function of the fraction p of connections, for three values of the 
membrane conductance over current ratio. The error roughly increases as ^/p, that 
is, the number of connections in the network. This scaling suggests that much of the 
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inference error is due to non-zero couplings. This finding agrees with Fig. [3p, which 
showed that the inferred interactions between uncoupled neurons was very small in the 
g — case. To better understand the performance of the algorithm, we compare in 
Fig. 2P the inferred interactions Jy with their true values for the 1560 oriented pairs 
j — s> i of a randomly drawn network of N = 40 neurons, with p — .2 and Jo = .2CVth ■ 
When the ratio gVth/I is small compared to unity, the quality of the inference is 
very good (Fig.|3p-1). For larger ratios gV t h/I the inferred couplings are still strongly 
correlated with their true values, but are approximately rescaled by an overall factor 
< 1, corresponding to the average slope of the linear regression in Fig. [4p-2. As gVth/I 
increases, this factor decreases and the inference error grows (Fig. |4j\) . 

Figure 3)3 shows that the inference error on the interactions increases not only with 
gVth/I but also with the noise ratio r. For large values of r, the network can sustain 
activity even when gVth > I, and the inference error can take large values (upper curve 
in Fig. [4)3). In this regime, the couplings found by the Fixed Threshold algorithm be- 
come small, and the inferred current Ii gets close to gVth- The corresponding potential 
V*(t) rises sharply, in a time r, to a value slightly below threshold, h/g, with small 
fluctuations due to the synaptic inputs. This phenomenon can be seen in Fig. |4p, 
which compares the optimal potential of a neuron for two different values of membrane 
conductance. As discussed in the Methods section, this behavior is a consequence of 
the a — > limit taken in the calculation of the optimal potential; when a, or r, is not 
small, the potential is unlikely to stay close to the threshold for a long time without 
producing a spike, see Fig. [2)3. In the next paragraph, we analyze the results of the 
Moving Threshold inference procedure. 

As a conclusion, zero couplings are perfectly inferred, while the amplitude of large 
(positive or negative) interactions can be underestimated by the Fixed Threshold al- 
gorithm, especially so when the noise is strong. However, the relative ordering of the 
interactions is essentially preserved by the inference procedure. 

3.1.4 Inference error with the Moving Threshold procedure 

The Moving Threshold procedure was tested in Fig. [2)3 on an asymmetric system of two 
IF neurons (J12/ (CV t h) = .1, J21 = 0) in the presence of a strong noise, see description 
in caption and Section 12. 61 While the Fixed Threshold procedure erroneously inferred 
that both interactions vanish, the Moving Threshold correctly inferred the sign and 
the order of the magnitude of the coupling: j ferred / '(CV ( /,) = .2 ± .1. The inferred 
currents were within 10% of their true values. These results were obtained from a large 
number S of spikes to avoid finite-S effects. 

The synthetic data used in Fig. [4)3 were generated with two different values of the 
noise ratio, r. We estimate the relative fluctuations of the potential around the optimal 
path, averaged over all the inter-spike intervals in the data set, using formula (| 14f) . and 
find 

y/{(Vj ~ gFj f -028 for r = .03 

V th 1 -138 for r = .15 1 ' 

for all values of gV t hl I comprised between .9 and 1.25. Hence, the relative fluctuations 
cannot be neglected when r — .15. Figure [4)3 shows the inference error obtained from 
the Moving Threshold algorithm as a function of the membrane conductance for that 
value of the noise ratio. Not surprisingly, the Moving Threshold procedure is more 
accurate than the Fixed Threshold algorithm. 
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In the Moving Threshold algorithm, the optimal potential is constrained to remain 
below a certain threshold, Vjjf , which depends on the time preceding the next spike 
and on the effective current If. Figure 0)3 shows the values of the moving threshold 
VM and of the optimal potential V* for a few spike intervals of the same neuron as 
in Fig. [4jD. As expected, the value of V*(t) lies substantially further away from the 
threshold V t h than in the Fixed Threshold procedure. In addition, Fig. [4f3 shows a 
random realization of the potential Vj(t), obtained through numerical integration of 
the LIF differential equation {TJ, for the same neuron i. Although Vi is stochastic, the 
comparison of several inter-spike intervals indicates that V*(t) and V,(i) are in fair 
statistical agreement. 

To investigate in more details the origin of the inference error on the couplings for 
large values of r and gVth/I, we plot in Fig. |4p the inferred values of the interaction 
Jij vs. the true value for every pairs j — > i of a randomly drawn network of N = 40 
neurons. The interactions inferred by the Fixed Threshold algorithm are about ten 
times smaller than their true values (Fig. [4p-3). The use of the Moving Threshold 
procedure leads to a spectacular improvement for positive- valued couplings (Fig.[4p-4). 
While positive couplings are accurately inferred, the magnitude of negative couplings 
is often overestimated. These negative couplings are responsible for most of the error 
e s in Fig. [4)3. From the Bayesian point of view, when r is smaller than the average ISI, 
negative-valued couplings are indeed intrinsically harder to infer than positive-valued 
ones. A positive input drives the potential closer to the threshold, which strongly 
reduces the ISI. Conversely, a negative input drives the potential down, and a spike is 
unlikely to occur before the potential first relaxes to its average level I jg after a time 
of the order of r. Hence, the influence of a negative input is hardly seen in the increase 
of the ISI when r is smaller than the average ISI. We present an analytical calculation 
supporting this argument in Section 2. 21 

3.2 Applications to multi-electrode recording data 

We now apply our algorithm to multi-electrode recordings of the ganglion cell activity 
of the salamander retina. Two data sets were considered. The first one, hereafter re- 
ferred to as Dark (data courtesy of M. Meister), reports the spontaneous activity of 32 
neurons for 2,000 seconds, and consists of 65, 525 spikes (Schnitzer and Meister, 2003). 
In the second experiment, referred to as Natural Movie (data courtesy of M. Berry), a 
retina was presented a 26.5 second-long movie, repeated 120 times, and the activity of 
40 neurons was registered for the whole duration of 3,180 seconds (Schneidman, Berry, 
Segev and Bialek, 2006). Natural Movie includes 172,521 spikes. The firing rates, av- 
eraged over the population of recorded neurons, have similar values in the two data 
sets: / ~ 1.02 spikes/sec in Dark, / ~ 1.35 spikes/sec in Natural Movie. 

These two data sets were analyzed in a previous work (Cocco, Leibler and Monas- 
son, 2009) with the perfect integrator model (g = 0) and the Fixed Threshold algo- 
rithm. In this section we extend the analysis to the case of the LIF model and use 
both the Fixed and Moving Threshold approaches. In particular we show that the LIF 
model is capable of inferring the asymmetry of the interactions, which is seen in the 
cross-correlograms but was not obtained with the perfect integrator model. Moreover 
we discuss error bars on the inferred couplings and the fact that strong negative in- 
teractions are more difficult to infer than positive-valued couplings. We stress that the 
couplings we infer a priori depend on the stimulus. Cocco, Leibler and Monasson (2009) 
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have studied how the interactions inferred with the perfect integrator model depended 
on the stimulus based on the analysis of two recordings on the same retina, namely 
the spontaneous activity and random flickering squares. An alternative approach to 
disentangle stimulus-induced and structural contributions to the couplings would be 
to consider a time- and stimulus-dependent external current I(t) fSection l4.4[) . 

The value of the membrane leaking time r strongly affects the number of contacts 
and the running time of the algorithm. It takes about 40 seconds to infer the currents 
and the interactions from either Dark or Natural Movie when r — 1 sec with one core 
of a 2.8 GHz Intel Core 2 Quad desktop computer, and about 10 times longer when 
t — 100 msec. The number of passive contacts of the optimal potential computed by 
the Fixed Threshold procedure quickly decreases as r increases. It is divided by ~ 20 
when the membrane leaking time increases from 100 msec to 10 sec for both data sets. 
In comparison, the number of active contacts is less sensitive to the value of r. We find 
that the ratio of the number of contacts per neuron and per second over the average 
firing rate takes similar values for both data sets. For r = 1 msec, this ratio is ~ 2.00 
for Dark, and ~ 2.04 for Natural Movie. The number of passive contacts is smaller with 
the Moving Threshold algorithm, while the number of active contacts remains rather 
unchanged compared to its value with the Fixed Threshold procedure. On the overall, 
the running time of the Moving Threshold procedure is higher due to the calculation 
of the time-dependent threshold . 

Knowledge of the variance of the noise is required for the Moving Threshold algo- 
rithm. The value of a could, in principle, be determined from experimental measures of 
the fluctuations of the synaptic current, but is unknown for the two recorded data sets 
available to us. We choose a so that the relative fluctuations of the potential around 
the optimal path V* are less than 10%. We compute these fluctuations by averaging 
(11411 over all ISI and all neurons i in the population. The corresponding value of the 
dimensionless standard deviation of the noise (|13p are: for Dark, a = .13, .12, .11 for, 
respectively, r = 200, 100, 20 msec; for Natural Movie, a = .15, .14, .12 for, respectively, 
r = 200, 100, 20 msec. 



3.2.1 Amplitudes of the inferred interactions and currents 



Figure [SJA. shows the average value of the current and of the interaction strength as a 
function of the membrane leaking time. As expected with the Fixed Threshold inference 
procedure, we find that the average value of the couplings decreases as r gets small. 
This effect varies from neuron to neuron: the closer ij is to gVth, the smaller are 
the couplings J.y . To compare the matrices of couplings J, J inferred with the Fixed 
Threshold algorithm for different values of r, we use the correlation coefficient (Hubert 
and Baker, 1979) 

f = cov(J.j') (26) 

Vcov(J, J) cov(J', J') 

where 

cov( J, f) = N(N - 1) J y - j'ij-iY, Jij) ( J2 4j) ■ (27) 

Identical matrices correspond to R = 1, and uncorrelated matrices give R — 0. R is 
independent of the scale of the coupling matrices J and J , i.e. R(aJ, a 1 J ) = R(J, J') 
for any a, a' > 0; therefore, R is sensitive to the relative amplitudes of the couplings 
J' and J and not to their absolute differences. We choose J to be the coupling matrix 
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Fig. 5 Amplitudes of the interactions and currents in Dark (full circles) and Natural Movie 
(empty diamonds). A. Average value of the current (left) and root mean square value of the 
coupling (right) as a function of the membrane leaking time r. Points corresponding to the 
Fixed Threshold (FT) procedure are joined by full lines, while dashed lines indicate the results 
from the Moving Threshold (MT) algorithm. Note that the currents are larger in Dark than in 
Natural Movie. B. Correlation coefficient R ill'lili between the couplings at leaking time r and 
with no leakage. C. Comparison between the interactions Jij found with the Moving (x-axis) 
and the Fixed (y-axis) Threshold procedures, for Dark (C-l) and Natural Movie (C-2). The 
dashed line is the x = y line, and r = 20 msec. D. Strongly negative couplings Jij vs. latency 
over the membrane leaking time r for three values of r. Couplings were obtained using the 
Moving Threshold procedure, and correspond to the Natural Movie data set. Only interactions 
Jij < —.1 are considered; there are, respectively, 16, 28, and 60 such couplings for r = 200, 
100, and 20 msec. The value of the slope of the best linear fit log(— Jij) = a latency(i, j)/t + /3, 
shown by the dashed lino, is a = 0.95. E. Distributions of the latencies (128ft between neurons 
in Dark (top) and Natural Movie (bottom). Only latencies larger than 5 msec are taken into 
account in the histograms. 



in the absence of leakage and J to be the coupling matrix for a given r. The value of 
R as a function of r is shown in Fig. [5)3. Even for r = 20 msec, the coupling matrix is 
substantially similar to the one obtained with the perfect integrator model (R — .6 for 
Dark, R = .5 for Natural Movie). Despite the overall change in the amplitude of the 
inferred couplings, the relative ordering of the couplings with the pair indices (i, j) is 
largely independent of r, especially so for Dark. However, for specific pairs of neurons, 
the interactions may strongly depend on r. Such a dependence effect will be illustrated 
in Section 13.2.31 and can be related to the temporal structure of the corresponding 
cross-correlograms . 

The average value of the interactions calculated by the Moving Threshold algo- 
rithm does not decrease when r gets smaller, and is larger than the one obtained from 
the Fixed Threshold procedure (Fig. [SJ\). To better understand this discrepancy, we 
compare in Fig. [5p the interactions inferred with both algorithms for every pairs of 
neurons in the Dark and Natural Movie data sets when r = 20 msec. The agreement be- 
tween both procedures is very good for positive and strong couplings. Couplings which 
are slightly positive with the Fixed Threshold procedure generally have a larger value 
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with the Moving Threshold procedure. This offset is responsible for the differences in 
the average values of the interactions found in Fig. [5)4.. In addition, in Natural Movie, 
negative- valued couplings often have a stronger amplitude with the Moving Threshold 
procedure. We find, in both approaches, a few negative and very strong couplings. The 
amplitude of those extreme couplings increases very quickly as the membrane leaking 
time decreases. 

The emergence of strong negative interactions with the lowering of r can be related 
to the presence of long latencies between the emission of spikes. We define the latency 
of neuron i with respect to neuron j as the smallest delay between a spike emitted by 
j and a later spike fired by i, 



A large value of the latency of neuron i with respect to j is interpreted by the inference 
procedure as the consequence of a strongly inhibiting coupling from j to i. However, the 
effect of a synaptic input of amplitude on the potential Vi of the neuron i decays 
exponentially with the ratio of the time elapsed from the input over the membrane 
leaking time. Hence, to keep the latency fixed while r is changed, the strong and 
negative interaction must change accordingly, 



where the constant has a negative value. Figure [Sp shows the negative couplings Jjj 
vs. the latencies of the corresponding pairs (i, j) divided by r, for three values of r. 
The outcome suggests that relation (|29|) is indeed correct, see Fig. [5p and its caption. 

The above mechanism explains why strongly negative couplings are less frequent in 
Dark than in Natural Movie. For r = 100 msec, there are 10 interactions (out of 1560) 
smaller than —1 in Natural Movie, and none (out of 992) in Dark. For r = 20 msec, 
these two numbers are equal to, respectively, 23 and 1. Figure [S}3 shows the histograms 
of latencies for both data sets. In Natural Movie, we find 17 pairs with latencies larger 
than 25 msec. In Dark, only one pair has a latency larger than 25 msec. The 

corresponding interaction, Jjj, is the only one smaller than —1 for r = 20 msec. 

3.2.2 Accuracy on the inferred interactions and currents 

As discussed in the Methods section, the uncertainty on the inferred parameters can 
be obtained from the Hessian matrix of L* , that is, from the curvature of the log- 
likelihood around its maximum. To quantify those uncertainties, we use the following 
procedure. Assume for instance we want to know how reliable is the inferred value, 
Ji,j , of the interaction Ji.j from neuron jo to neuron i. We fix Jij to an arbitrary 
value, and maximize L*(T\ Jij, 1%) (|16p over all the couplings with j ^ jo and over 
the current Ii . The outcome is a function of Jij , which we denote by L c and call 
marginal log-likelihood. L c (Jij ) has, by definition, a maximum in Jij = Jij - Its 
second derivative in the maximum, L"(^i,i )' * s related to the error bar AJi j on the 
interaction through, see (|17|) and l|18|). 



latency(i,j) 



mm 

k,£:t i-k <tje<ti tk + 1 



(*i,fe+l — tj,i) 



(28) 




(29) 




(30) 
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Fig. 6 A. Marginal log-likelihoods L c (Ii) (top panel), and Z/ C (^i,27)> L c (Ji 20)1 L c (Jii) 
(from left to right in the bottom panel) for Natural Movie, and r = 200 msec. Dashed lines 
correspond to the best fits with a single quadratic function. The most likely value for the current 
is I\ = 1.14gl/ t h. The most likely values for the interactions arc: Ji,27 = —.11, Jl,20 = -01, and 
Ji,4 = .22, in units of CV t h- The offset on the vertical axis has been chosen so that all maxima 
are at height L c = 0. B. Average value of the first-passage time tfPT after a synaptic entrance 
of amplitude J. C. Derivative of tfPT with respect to J. The parameters of the neuron are: 
gVth/I = 1-2, r = .15, r = 85 msec (full line) and 20 msec (dashed line). The derivative is 
maximal around Jopt/{CVth) = 1 — I/(gVth) — -167. 



The same procedure can obviously be used to obtain the error bar on the current 7j. 

We now illustrate this approach on the Natural Movie data set, and one arbitrarily 
chosen neuron, i = 1. Three interactions, representative of, respectively, positive, weak, 
and negative couplings, were singled out among the 39 couplings incoming onto neu- 
ron 1. Figure[6j\ shows the marginal log-likelihoods L C {J\^), ic(Ji,2o)i an d L c (Ji, 27)1 
in addition to L c (Ii). For all four parameters, the marginal likelihoods can be approxi- 
mated with parabolas in the vicinity of their maxima. Estimating the second derivatives 
from those best quadratic fits and using (|30p . we obtain 

Ah AJi 27 AJi 20 AJi 4 . , 

where a is the dimensionless noise level defined in (|13|l . Hence, the error bars on the 
couplings and currents have very similar values. This common value depends on the 
noise level, a. As discussed in the next Section [3.2.11 a is expected to be close to, or 
smaller than unity when r = 200 msec. Consequently, the value for Ji,20 is compatible 
with zero, while the interactions J\ 27 and are non zero, with 99.9999% confidence. 

A closer inspection of Fig. [SJA. shows that the quality of the quadratic fit of L c is 
excellent for Ji 27 an d J\ 4, but less so for I\ and Ji,20- For the latter parameters, 
it seems that the curvature of L c takes two different values, depending on whether 
the maximum is approached from the left of from the right. This phenomenon results 
from the piece-wise structure of the L* function, see Methods section. A practical 
consequence is that the errors I\ — I* and Ji^o — jl,20 are n °t evenly distributed 
around zero; for instance Ji^o is more likely to be larger than Ji^o than it is to be 
smaller. 

Note that strong, negative interactions may be harder to infer than positive- valued 
couplings, a phenomenon already underlined by Aersten and Gerstein (1985). The 
underlying intuition is that the duration of the ISI is less affected by an inhibitory input 
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Fig. 7 A. Cross-corrclograms H(t) for pairs (5, 17) and (1, 22) in Dark. The cross-correlograms 
are normalized such that H(t) — > 1 for large delays \t\. B. Ratios Jij/Jji of the interactions 
between the neurons 5,17 (top) and 1, 22 (bottom) as a function of r. 



than by an excitatory input when the membrane leaking time, r, is small compared to 
the average value of the ISI. We now present an analytical argument supporting this 
intuition. Consider a neuron, fed with an external current / and with noise variance 
equal to a 1 . Assume a synaptic input of amplitude J is received at time t = 0. We 
call tppx the average value of the time at which the neuron will emit a spike; the 
calculation of tppp can be done using a series of parabolic cylinder functions (Alili, 
Patie and Pedersen 2005) . Figures |6j3&C shows that the dependence of tp pp on J is 
much weaker for negative- valued J than for positive couplings. As the set of spiking 
times is the only information we have at our disposal, the difficulty in inferring negative 
couplings is intrinsic to the Bayesian approach, and cannot be circumvented by any 
particular algorithm. 

3.2.3 Symmetry of the interactions and cross-correlograms 

The dependence of the symmetry of couplings upon the membrane leaking time r can 
be understood, to some extent, from the structure of the cross-correlograms, that is, 
the histograms Hij(t) of the delays t = — tjj between the times of the spikes fired 
by the two neurons i, j in each pair. To do so, we consider two pairs of neurons in 
Dark, called pairs (5,17) and (1,22). Figure [7J\ shows the cross-correlograms H^^yj 
and H\ 22- Pair (5, 17) is characterized by a positive peak in H , centered in t = 0, and 
of width ~ 20 msec. Pair (1,22) exhibits a positive peak of correlations, of the same 
width, but centered around t ~ 20 msec. 

We plot in Fig. [7)3 the symmetry ratios of the interactions in the pairs, P5,i7 = 
•^5,17/^17,5 and Pi,22 = ^1,22/^22,1 • We find that £5,17 is, to a large extent, independent 
of t. Conversely, pi.22 sharply decreases with decreasing r and is close to zero when 
t — 20 msec, which coincides with the typical delay in the cross-correlogram H\ 22 
shown in Fig. [7J\. We conclude that the inference procedure is capable of capturing 
the directionality of the interaction between the neurons 1 and 22, if r is small enough. 
This results shed some light on the correspondence between the interactions inferred 
within the LIF model and within the Ising model (Schneidman, Berry, Segev and 
Bialek, 2006; Shlens et al, 2006). Couplings inferred with the perfect integrator model 
for Dark are in good agreement with the Ising interactions, when the time is binned into 
windows of width At = 20 msec (Cocco, Leibler and Monasson, 2009). By construction, 
the Ising model produces symmetric interactions from the pair-wise correlations of the 
activities, averaged of the binning window. In the absence of leakage, the Integrate-and- 
Fire inference algorithm hardly distinguishes between a post-synaptic and pre-synaptic 
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firing pattern, and produces rather symmetric couplings. But as r decreases, the LIF 
couplings may become strongly asymmetric (Fig. [7)3). In this case, the correspondence 
between the Ising and LIF couplings breaks down. The same phenomenon was observed 
in Natural Movie, where delays in the cross-correlograms are even stronger. 



4 Discussion 

In this article, we have presented a procedure to infer the interactions and currents 
in a network of Leaky Integrate-and-Fire neurons from their spiking activity. The 
validity of the procedure was established through numerical tests on synthetic data 
generated from networks with known couplings. We have also applied our algorithm 
to real recordings of the activity of tens of ganglion neurons in the salamander retina. 
Though our algorithm is limited to moderate noise levels and instantaneous synaptic 
integration, it is fast and can, to our knowledge, handle much bigger data sets than 
the existing inference methods for the stochastic IF model. It is our intention to make 
this algorithm available to the neurobiology community in a near future. 



4.1 Comparison with previous studies 

Cross-correlation analysis (Perkel, Gerstein and Moore, 1967; Aersten and Gerstein, 
1985) consists in studying the distribution of delays between the spikes of neurons in 
a pair. This approach has been used to characterize the connections between neurons 
(amplitude, time-scale, dependence on distance), or their dynamical evolution (Fuji- 
sawa, Amarasingham, Harrison and Buzsaki, 2008). The analysis do not require any 
combinatorial processing of the activity of a large part of the neural assembly. As a re- 
sult, the approach is not limited to small networks. However, cross-correlation analysis 
may find difficult to separate direct correlations from indirect correlations modulated 
through interactions with neurons in the surrounding network (Ostojic, Brunei and 
Hakim, 2009; Cocco, Leibler and Monasson, 2009), or due to common inputs (Con- 
stantidinidis, Franowicz and Goldman-Rakic, 2001; Trong and Rieke, 2008). 

In statistical approaches a widely-used concept is the one of causality (Seth and 
Edelman, 2007). A causal interaction exists from neuron i to neuron j if the knowledge 
of the activity of i helps predict the firing of j beyond what can be achieved from 
the activity of j alone. In practice, causal relationships are detected through linear 
multivariate statistical regressions (Sameshima and Baccala, 1999), and may overlook 
non-linear dependencies. Causal analysis have also difficulties in evaluating the strength 
of the interactions. 

Maximum entropy models, which deduce interactions from pairwise correlations 
only, have been shown to accurately reproduce higher-order correlations between neu- 
rons in the vertebrate retina (Schneidman, Berry, Segev and Bialek, 2006; Shlens et al, 
2006; Cocco, Leibler and Monasson, 2009). These models, however, suffer from some 
limitations. Interactions are constrained to be symmetric, and temporal correlations 
are partially discarded (Marre, El Boustani, Fregnac and Destexhe, 2009). In addition 
obtaining the interactions from the correlations may be computationally very hard for 
large networks, though efficient approximate algorithms have recently been developed 
(Cocco and Monasson, 2010). 



21 



Generalized linear models (GLM), which represents the generation of spikes as 
a Poisson process with a time-dependent rate, have been applied to various neural 
systems (Brown, Nguyen, Frank, Wilson and Solo, 2001; Truccolo et al, 2005; Pillow 
et al, 2008). The inference of parameters in the GLM framework is apparently easier 
to solve than for IF models, which has made the GLM framework very attractive. 
Whether GLM are better than IF models to account for real neural activity, regardless 
of the computational complexity of both inference framework, is an important issue 
(Gertsner and Naud, 2009). We hope that our work, which makes possible to apply the 
IF model to large data sets, will help to answer this question. 

Approaches to infer model parameters in the IF framework have been so far capable 
of processing a very limited number of neurons or of spikes. Pillow et al. (2005) inferred 
the model parameters of one stochastic IF neuron based on a 50 second-long recording 
with a procedure tolerating any level of noise; Makarov, Panetsos and de Feo (2005) 
inferred the connections between 5 deterministic IF neurons from a 60 second-long 
synthetic spike train. In comparison we have analyzed a 3180-second long recording of 
the activity of 40 neurons. 

The running time of our procedure increases as N 2 S ~ N 3 T /, where T is the du- 
ration of the recording and / is the average firing rate. Recently, Koyama and Paninski 
(2009) have proposed a numerical procedure for calculating the optimal potential and 
inferring the interactions. In their approach, the time is discretized into many time- 
bins of small duration A, and the values of the optimal potentials at those discrete 
times can be found by means of the interior-point method for discrete constrained opti- 
mization problems. The running time of the procedure, 0(N 3 T/A), is approximately 
l/(fA) times larger than ours. In practice, / is of the order of 1 to 10 Hz, while the 
discretization time, A, is of the order of 1 msec; hence, l/(fA) ranges from 100 to 1000. 
However, this order of magnitude does not take into account the existence of multi- 
plicative constants; a comparative test of the two approaches on the same synthetic or 
real data would be useful to accurately estimate their running times. Furthermore, the 
algorithm introduced by Koyama and Paninski can easily incorporate the presence of 
temporal filtering in the interactions. Our procedure is, in its present form, valid when 
the integration kernel is instantaneous only; considering other synaptic kernels would 
require ad hoc modifications to the expressions of the optimal noise and potential and 
to the search procedure for contacts. 

4.2 How to include a finite integration time 

One of the major assumptions in our approach is that the synaptic integration time, 
r s , is vanishingly small. In practice, r s does not vanish, but might often be smaller 
than the membrane leaking time, r, and the average ISI. Assume that neuron i, whose 
potential Vi is close to the threshold Vj^, receives a spike at time t from another neuron, 
j, through a strongly excitatory connection J,j > Vth — Vi. Then, neuron i will reach 
the threshold level after having received a charge Aq — C(Vth — Vi), smaller than Jy. 
As a consequence, large positive interactions can be underestimated when the latency 
of neuron i from neuron j (|28l) is smaller than r s . 

To compensate for this effect we could introduce a time-dependent value for the 
interaction, 




(32) 
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where tik+i(> t) is the closest firing time of neuron i. Hence the effective interaction 
Jij{t,~ti,k+l) is equal to its nominal value Jij only if the synaptic current has enough 
time to enter the neuron j, and is a fraction of otherwise. The modified procedure 
will be correct as long as r s < r. If the synaptic and membrane time-scales are com- 
parable, one needs to take into account the complete shape of the synaptic integration 
kernel, K(t). Choosing simple enough integration synaptic kernel, such as the piece- 
wise linear function K(t) = if t < or t > t s , K{t) = 2 mm(t, t s - i)/rf if < t < r s , 
could lead to tractable dynamical equations for the optimal potential and noise. The 
resolution of those equations is left for future work. 

4.3 Towards a more realistic inference model 

The inference procedure that we have introduced here can be extended to include 
realistic features such as a refractory period, tr. To do so, we restrict the sum in (11011 
to the spikes m entering the neuron i at times larger than to + T R- We have run the 
modified inference procedure on the recordings of the retinal activity, for values of tr 
ranging from 2 to 5 milliseconds. The couplings did not change much with respect to 
the values found with tr = 0. Note that the introduction of a propagation delay tjj in 
the synaptic interaction is straightforward, as long as the integration kernel remains a 
Dirac distribution (centered in rrj). 

Bounds on the values of the couplings and currents e.g. to prevent the exponen- 
tial growth of negative interactions with the leaking conductance can naturally be 
introduced through a prior distribution. As an example, assume that the interactions 
Jij take values in [J_, J+]. Then, one could maximize L* — ^""^ W( Jjj) instead of the 

i-.j 

log-likelihood L* alone, where W(J) = f (J - J_ f if J < J- , if J_ < J < J+, 
f (J - J+) 2 if J > J+ and to is a large positive coefficient. 

We have assumed, throughout this work, that the values of g and Vj^ were known. 
In practical situations, while the orders of magnitudes are known, the precise values of 
these parameters should be inferred, and could depend on the neuron i. The inference 
procedure could be modified to update the values of gi and (Vth)i a t f ne same time 
as the synaptic couplings Jij and the current 7^. The number of parameters to infer 
(per neuron) would simply increase from JV to iV + 2, and the running time should not 
increase too much. 

4.4 Inference from a limited neural population and in the presence of a stimulus 

Nowadays, multi-electrode experiments can record a few tens, or hundreds of neurons. 
To which extent do the interactions inferred from this sub-population coincide with 
the interactions one would find from the knowledge of the whole population activity? 
The question does not arise in cross-correlation analysis: the correlation between the 
firing activities of two neurons is obviously independent of whether a third neuron is 
recorded or not. However the issue must be addressed as soon as a collective model for 
generating the activity is assumed, such as the coupled LIF models studied here. 

A detailed analysis suggests that the interaction between a pair of neurons is not 
affected by the activity of other neurons distant by more than I = 300 \xm in the 
case of spontaneous activity (Cocco, Leibler and Monasson, 2009). The electrode array 
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should be at least twice longer and wider than £, and should be dense enough to 
capture all the neurons on the recorded area. It is estimated that about 10% of the 
ganglion cells are registered in the Dark experiment, compared to more than 80% with 
the denser but smaller electrode array used in the Natural Movie experiment (Segev, 
Puchalla and Berry, 2005). It would thus be very interesting to repeat our study on 
other multi-electrode recordings, with sufficiently large and dense arrays. 

Taking into account the stimulus S in the inference process would also be inter- 
esting. To do so, we could add a stimulus-induced current, I s (t\S), to ([T}. A simple 
expression for this current would be J s (f|S) = J dt'Kf(t — t) Si(t'), where K" is a 
kernel similar to the one used in generalized linear models (Pillow et al, 2008) . The ex- 
pression of the current-dependent term in the potential V(rj, t) (|10|) should be modified 
accordingly, while the noise-dependent term would remain unchanged. It is important 
to note that the search procedure for contacts presented in Section [2.41 would remain 
valid. However, the expressions of the noise coefficient, the contact time and the du- 
ration of a passive contact given in Appendix [B] for the case of a constant current 
/ should be rederived and would depend on the precise temporal structure of the 
stimulus-induced current I s (t\S). 
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A Active contacts 

In this Appendix, wc justify the prescriptions in the search for active contacts presented in 
Section 12.41 For the sake of simplicity we restrict to the g = case (no membrane leakage); 
the extension to non-zero g is briefly discussed in Appendix [B] We consider a neuron i, and 
call M the number of spikes received by this neuron during its k th inter-spike interval [to = 
ti k \ *M+l = *i fe+i]- The arrival times are t% < t% < . . . < tjvfi an d the corresponding synaptic 
strengths are Ji, Js, . . . , Jm- To lighten notations we hereafter omit the index i of the neuron. 

A.l Case of M = or 1 input spike 

To understand the key notion of contact, we first consider the simple case of a neuron receiving 
no spike during the inter-spike interval [tj j. = 0; ti, fc+i = T]. The optimal noise is constant 
according to 10. Equation (0 then shows that the optimal potential is a linear function of the 
time, which is fully determined from the boundary conditions V*(0) = 0, V*(T) = Vth- We 
obtain 

V*(t)=V th ± and „«( t ) = £p* - J . (33) 

This solution is correct since the potential remains below the threshold at all times < t < T. 

Let us now assume now that the neuron receives one input from another neuron, of strength 
Jl at time t\ £]0; T[. The effect of the input is a discontinuous jump of the potential at time 
t\ and of size -jj, shown in Fig. [H] Repeating the calculation above, we obtain the following 
expressions for the optimal potential and noise 

VX(t)=(v th -£\±+^e(t-tt) and V * A = CVth ~ Jl - I (case A), (34) 
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Fig. 8 Sketches of the optimal potentials V* (top) and noises 77* (bottom) for one neuron 
receiving one weak (A), one strong negative (B), and one strong positive (C) input. The 
jump in the optimal noise consecutive to an active contact is always positive. Values of the 
parameters used for the figure are: J = 0, ti = T/2, J 1 /(CV th ) = .2 (A), -1.2 (B), 1.2 (C). 



where 9 is the Hcavisidc function: 9(x) = 1 if x > 0, otherwise. This solution is sketched in 
Fig. [S^. It is valid when the potential Vjf remains below the threshold at all times. We call 
this situation case A. As V A is a piece- wise linear function we only need to check that V^J(tj~) 
and Va(*i~) are both smaller than Vth- The two conditions are fulfilled provided that 



J- 



-CV, 



T-ti 
1 ti 



<J!<J+ = CV t , 



(35) 



What happens when the above condition is violated? Let us consider first J\ < J_ (referred 
to as case B hereafter). Then Vjf exceeds the threshold Vth before the input enters the neuron. 
To prevent the potential from crossing the threshold at time ti, the true optimal noise, T]g, 
should be smaller than r) A . But, if r)* B < r)* A , the potential could not reach Vth when the neuron 
emits its spike at time T according to the very definition of rj A \ The only way out is that r] B 
takes two different values corresponding to the two sub-intervals [0;ti[ and ]ti;T], which we 
call, respectively, if B _ and tj b We expect rj* B _ < rj* A < rf g ,. The noise can change value 
in i = ti through J9} only if the potential reaches the threshold in ti . We find that 



vS(*) = v« 



'1 



and 



cv tl 
tl 



(case B, < t < ti) 



(36) 



from the boundary conditions V*(0) = 0, V*(t 1 ) = Vth, an< i 
Ji T-t , , Ji 



V B (t) = V th + 



C T-t! 



and 



Vb, 



(case B,ii < t < T) , (37) 



from the boundary conditions V*(t~^) = Vth + ~Q~tV*(T) = Vth- This solution is drawn in 
Fig. [HP- It is important to stress that the above solution is based on the capability of the noise 
to abruptly change its value when the potential touches the threshold in t = ti. A detailed 
study of the behavior of the noise close to such 'contact points' proving that this is indeed the 
case is postponed to Appcndix lA.2l 

Finally, we turn to case C corresponding to Ji > . In this case the input is so excitatory 
that the noise has to be negative to prevent the neuron from emitting a spike at a time t < ti . 
As in case B, the potential reaches the threshold in t = ti to allow the noise to change its 
value after the input has entered the neuron. We find 



vS(t) 



Vtfc - § 1 



and 



CVth - Ji 

tl 



(case C,0 < t < ti) , (38) 



according to the boundary conditions Vq (0) = 0, Vq (t-^ ) = Vth — q ■ Right after the spike has 
been received, the potential has reached its threshold value, and will keep to this value until a 
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spike is emitted at time T, hence 

Vg{t) = V th and Tfc + = -I (case C, ii < t < T) . (39) 

This solution is drawn in Fig.[8p. 

We now give the values of log-likelihoods L* corresponding to the cases listed above. The 
value of L* can be calculated from the knowledge of the optimal noise 77* through JTBJ). In the 
case of M = spike, we find, using 11331 1 with T = t\ — to, 

L {UhtllI) = — ■ (40) 

The optimal current is then inferred by maximizing L*(I) with the result I = tl L to > which 
corresponds to a vanishing value for the optimal noise, as expected. 

When M = 1 spike is received by the neuron, the log-likelihood L* has three distinct 
expressions corresponding to the case A, B, C discussed in Section fA.il The resulting expression 
is (with t 2 =T): 

- (CVt ""2fe"- t or to))2 if J-<Jl< J+ (^se A) 

- (CVth ~2fa~- t ( or to))2 - - fc) if Ji > j + ( case °) 

(41) 

The log-likelihood L* is a continuous and convex function of its argument. The first derivatives 
of L* arc continuous in J_, J+, but the second derivatives are not. 



A. 2 Study of the optimal noise close to an active contact point 

The noise coefficient n in ifSJl arc constant over the time interval separating two active contacts. 
The value of r\ may however change upon the crossing of an active contact. The scope of this 
section is to show that the noise right after the contact can take any value larger than the 
noise immediately before the contact. This monotonicity property justifies the search for the 
minimal noise coefficient done in (lilt , see Appendix IA. 31 

To show that the noise always increases through an active contact, we consider that the 
synaptic integration is not instantaneous, but takes place over a finite albeit small time, t s . 
We thus replace the expression for the current J?" in ||2} with 

where J;j is the strength of the connection from neuron j onto neuron i, and K(t) is is the 
memory kernel of the integration of synaptic entries (top panel in Fig. [9j . We assume that 
K(t) vanishes for r < and for r > t s where the integration time t s is independent of the 
pair In addition, K is positive, and its integral over the interval [0;r s ] is equal to unity. 

We consider the case of a single incoming spike, as in Scction lA.il We want to show that, 
in the t s — > limit, the only constraint linking the values 77* and 77^ of the optimal noise, 
respectively, before and after a spike entering at ti , is rf^ > 77* , as we have found for a single 
incoming input in cases B and C. To do so, we assume that the time of synaptic integration r s 
is small but finite , and consider case B. The dynamics of V* and 77* can be divided in several 
steps, whose numbered are reported on Fig. [9] 

1. Prior to the input, i.e. at times < ti, the optimal noise 77* is constant and the optimal 
potential V* is a linear function of the time, with slope {I +n*_)/C , as shown in Fig.(9fleft). 

2. A strongly negative input of amplitude Ji(< J-) is then received by the neuron between 
times ii and t\ + r s . The derivative of the potential now decreases with the time until it 
vanishes at time t c defined through 

n* +1 

K(t c -ti) = k- where k- = . (43) 

—Jl 
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Fig. 9 Behaviors of the optimal potential V* (middle) and noise 77* (bottom) close to a 
contact point, compared to the memory kernel K (top). An input of total amplitude Ji enters 
the neuron during the time interval t\ < t < ii +r s . Left: J\ is strongly negative as in Fig.[8jS; 
italic numbers refer to the steps listed in the main text. Right: Ji is strongly positive as in 
Fig.[5p. Sec text for a detailed description of the curves, of the constants (1431441) . and 

of the times t\,t c ,t' c , f", t s . 



3. If the value of 77*_ is chosen so that V*(t c ) = V t h, the potential tangentially reaches the 
threshold at t c (contact point). Then, the potential remains constant and equal to Vth- 
The noise obeys eqn. ||9} and, therefore, increases until it reaches the prescribed value, rA, 
at time t' c such that 

vl+i 

K(t' c -t 1 ) = k+ where k+ = -!± , (44) 

— Jl 

sec bottom panel in Fig- [9lleft). 

4. Then the potential starts decreasing from its threshold value through eqn. (|6j, and reaches 
a minimum in t", solution of the same equation i !44[ ) as t' c , see Fig. [9f top left). 

5. At later times the derivative of the potential is positive from eqn. JB), and increases until 
time ti + t b , coinciding with the end of the synaptic integration. 

6. At times larger than ti + t s , the potential keeps growing with a constant slope equal to 

(i + vl)/c 

In the T a — > limit, all times t c ,t£.,t" tend to the same value, that is, the time ti. More 
precisely, as the slope of K is of the order of t s ~ 2 (in absolute value), and 77* , 77^, V* (ii) are 
finite (= O(l)), then for t s — > 0, ti,t c ,t' c differ from each other by O(r^). Hence the change 
in the potential V* between t' c and t\ + t s equals -^4- + 0(t s ). We conclude that, for t s — > 
the potential becomes a discontinuous function of time with a discontinuity ^r. In addition, 
the noise rf can also jump abruptly from its value rf_ at to any larger value 7/1 at time 
since the maximum of K tends to infinity when t s — > 0. 

Note that the drawing of Fig.^left) tacitly assumes that fc+ > A hypothetic scenario 
would be that the noise exactly compensates the synaptic input for a longer time interval 
(including the top of K), while the potential would remain equal to Vth - I n this case, the peak 
value of the noise would be 0(l/r 3 ). The contribution to the integral (1161) would be of the 
order of 1/t s and would diverge in the t s — > limit. Hence this possibility is precluded. 



30 



The above discussion is straightforwardly extended to case C. The optimal potential and 
noise are sketched in Fig. [9f right). Note that the contact interval spreads beyond [t' c ,t c ] in this 
case. In the generic case of more than one incoming spikes, the contact interval is restricted to 
[tci^c] as in case B. The noise can also discontinuously change from its value rf_ < —I before 
the contact to any larger value, rf, , after the contact. 



A. 3 Case of M > 2 incoming spikes 

We now consider the general case of M input spikes. Let Vq = 0, mo = 1 be, respectively, the 
initial value of the potential and the index of the first input spike. We define the piece-wise 
linear function solution of {6]l for a constant noise rj, 

j , M j 

V(r,,t,t ) = V + ^p-(t-to)+ ]T J^e{t-t m ). (45) 

m— mg 

We are looking for the smallest value of the noise coefficient r\ capable of bringing the potential 
V(i), t, to) from its initial value V(rj, to, to) = to the threshold. The contact time, t c , coincides 
with an entering spike, i.e. t c = t m * for some m* > 1. If m* = M+l then the optimal potential 
is V(rj* ,t,to) throughout the inter-spike interval [to;tM+i], and the problem is solved. If 
m* < M, t m * is the first active contact point of the potential. rf and V(rj* , t) are, respectively, 
the optimal noise and potential on the interval [to,t m *]. 

The correctness of the above statement can be established using a proof by contradiction. 

— assume that the optimal noise, r] opt , is smaller than rf on some sub-interval of [to;t c ]- 
Remark that the potential V in (1456 is an increasing function of the noise, 

rf >r?=> V(r/,Mo) > V{r,,t,t ) , (46) 

for all t > to- By virtue of 1461 1 and the definition of rj* , V(r]° pt , t, to) cannot touch the 
threshold at any time so the noise is constant throughout the interval [to;t c ]. Hence no 
active contact can take place at time t c . As rj* is the minimal value of the noise which 
can drive the potential into contact with the threshold over [io;*M+lL we conclude that 
V{rj opt ,t, to) cannot cross the threshold at any time < <m+i- The neuron can therefore 
not spike at time tM+i- 

— conversely, suppose that the optimal noise is equal to r] a > 77* on the interval [foitm ] 
with 1 < m a < m* , and takes another value on the interval [t m a • t c Q. As the change of 
noise can take place only through an active contact, and the change is necessarily positive 
(Section I A. 21 . we have ry^ > rf . Applying 1146 II to the interval [t m a- t t c ], we have 

V {rf ,t c ,t m c)> V( V a , t c , t m c ) . (47) 

Adding the value of the optimal potential in t m a to both members of the previous in- 
equality, we find 

V*(t c ) = V(r 1 f 3 ,t c ,t rn «) + V{ V a ,t nia ,to) 

> V(r 1 a ,t c ,t ma ) + V(ri a ,t rn «,to) 
= V{r, a ,t c ,t ) 

> V(ri*,t c ,t ) (48) 

where the last inequality comes from (146 1 ■ But, by definition of rj* , V(rf ,t c ,to) = Vth- 
Hence, we find that the optimal potential in t c is above threshold, which cannot be true. 
The optimal noise and potential on the remaining part [t c ; £m+i] of the inter-spike interval 
can be determined iteratively. We replace to with i m » and Vo with V t h if Jm* > or Vth + 
if J m * < in 1451 . and look for the lowest noise producing a new contact point over the 
interval [t m » , tftf_|_i]. The procedure is repeated until the whole interval is exhausted. This 
way an increasing sequence of noise values is obtained, each corresponding to the slope of the 
optimal potential between two successive contact points. 



7 The case of three or a higher number of values for the noise can be handled exactly in the 
same way. 
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Fig. 10 Sketch of the optimal potential close to a passive contact starting at time t c . The 
duration of the passive contact is A c (£), where t is the index of time tg corresponding to the 
next active contact. The potentials corresponding to the hypothesis 1=1 and 1 = 2 are shown 
with the dashed and full curves respectively. 



B Passive contacts 

When the membrane leaking conductance is non zero, some change have to be brought to the 
above calculation of the optimal noise and potential. First, in the absence of inputs, the noise 
is no longer constant, but rather it is an exponentially increasing (in absolute value) function 
of the time Similarly, the potential V itself is not a linear function of the time as in II 181) . 
but is a linear combination of cxp(±t/r) with appropriate coefficients, sec HOI) . 

The main conclusion of AppcndixfX]still holds: the difference between the noise values just 
after and before an active contact point, coinciding with a synaptic input, is always positive 
(Fig. [2}\). Consequently, the procedure of Section 12.41 i.e. the iterative search for the active 
contact points and the minimal noise coefficient rf , defined through Jill) , remains unchanged. 
Note that some care must be taken to translate the statement about the growth of the noise 
to the values of the noise coefficients. Consider for instance two successive contact times, t and 
t' , and call r\, r/ the corresponding noise coefficients. That the noise is larger at time t' than 
at time t implies that 77 X exp((4' — t)/i~) < rf , but does not imply that r/' is larger than r/[f]. 

There exists, however, a major difference between the g = and g ^ cases. When 
g > 0, the optimal potential is not guaranteed to be a monotonous function of the time, as 
shown in Fig. 1101 For given values of g, I, and of the times and the amplitudes of the synaptic 
inputs, the optimal potential V* may touch the threshold at an intermediate time, t c . Such a 
situation is called passive contact. It is important to note that the value of the optimal noise 
during a passive contact is, according to eqn. equal to gV t h ~ !• As the optimal noise is a 
monotonous function of the time between two active contacts, see eqn {HJ, the value gVth — I 
can be crossed at most once: there is at most one passive contact in between two successive 
active ones. To be more precise, there are at most A + 1 passive contacts in an inter-spike 
interval with A active contacts. 

To decide the existence of a passive contact in an interval [to; t.M+l]> we look for a solution 
of the two coupled equations expressing that the optimal potential touches the threshold 
without crossing it, 

8V* 

V*(v P ,t c ) = V th and -— (?7,tc) = . (49) 
at 

The solutions of these equations give the noise coefficient r\ p and the contact time t c at which 
the optimal potential reaches the threshold value (Fig. HUH . The solution can be calculated 
analytically, with the following result. Let us call Vo the value of the potential of the neuron 
at time tg . For each m < M we define 

V m = V + J2 ^e^-*o)/r t (so) 

t:t <t t <t m 



8 This situation can not happen in the 9 = case, where the noise and the noise coefficient 
coincide. 
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where the summation runs overs the spikes entering the neuron between times to and t c . A 
passive contact takes place in the interval [t m ;4 m _|_i] if: 



gVth ~ I an d Vm — Vth have the same sign; 
the noise coefficient 



ri P = gV m - I - ^/{gVm - I) 2 - (gV th - If 



(51) 



is smaller than the lowest noise coefficient corresponding to all the possible active contacts 
at times t e , with 1 < £ < M; 
the corresponding contact time 



to — t log 



Vp 



gv th 



(52) 



where rj p is given by ( 1511 . lies in the correct interval: t m < t c < £ m +i; 
the optimal potential can reach again the threshold at a later time, coinciding with an input 
spike or with the end of the inter-spike interval. We call A c (£) the duration of the active 
contact such that the potential reaches V t h a t time ti, starting from Vth a t time t c + A c (l), 
see Fig. 1101 The analytical expression for the duration of the passive contact allowing the 
potential to be in active contact at time tt is 



= -TlOg 



2V a (£) 



V b (i) ■ 



(53) 



whore 



V a {i) = ^ e ^"*o)/ T and V b {l) = V th 

2g 



g , c 



Jt 
c 



9(Jt) ■ (54) 



We must have t c + A c (i) < ti for at least one value of i > m + 1 . 

When all the conditions are fulfilled, a passive contact is present. The duration of the contact, 
A c , merely plays the role of a latency time after which the potential V* resumes its course 
(Fig llOII . We can check that V* is an increasing function of A c . The optimal value of A c will 
therefore be equal to the smallest possible value of A c (£), for the very same reason that we 
had to chose the minimal noise when looking for active contacts, see example in Fig. 1101 

To end this Appendix, we give the expression for the log-likelihood L* 11161 for an interval 
including a passive contact between two active contacts. Gathering the contributions to the 
integral of the squared optimal noise coming from the three intervals [to; t c ], [t c , t c + A c ], and 
[t c + A c ; t m * ] , we obtain 



L*(T\J,X) : 



(gV th - If 



exp [2(t m , -t c - A c )/rj - exp [ - 2(t c - t )/r] 



(55) 

Differentiation of L* with respect to the current and couplings gives the expressions for the 
gradient and Hessian matrix needed for the Newton-Raphson method. The expressions are 
easy to obtain but are lengthy, and thus we do not reproduce them. 



C On the eigenmodes of the Hessian matrix for weak couplings 

In this Appendix, we analyze the eigenmodes and eigenvalues of the Hessian matrix of the 
log-likelihood L*, and relate the eigenmodes to the fluctuations of the effective current, If, of 
the current, Ii, and of the couplings, Jij. Consider two successive spikes emitted by neuron 
i and the optimal potential V*(t) on the time interval [tj fc;tj fc+i]- When the couplings Jij 
vanish and passive contacts are absent, V*(t) does not enter into contact with the threshold 
at times < tj fc+i. By continuity, this statement remains true if the couplings Jij have very 
small amplitudes. In this regime, the stochastic process undergone by the potential is simply 
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the Ornstein-Uhlcnbcck process with a time-varying force, and the expression for L* (116ft is 
exactly given by 



L*{T\J,1) = - 1 -J2^ ) ( CV ^~ E J i3^-hr<i>%) (56) 

i,h V j&i) 



where 



and 



(0 = 1 (i_ e -2(ti, fc+ i-ti, fc )/r) t (57) 



>8 



(58) 



; 

1 _ e -(*i,fc+i-*i,fc)/T ifj = i. 

The Hessian matrix of L*, attached to neuron i, is the N X N matrix I I17II with elements 



HW is a positive matrix according to J 59 1) - To study its spectrum let us first consider the case 
)ry i 

the k th ISI of neuron i, we have 



of very weak leakage (very large t). In this limit, calling At£ = tj^+i — ti,jfc the duration of 



(i) 1 fi) ^ifc (i) ih 

~ * (7y ' ^k i ~ * ' ^fc j ~~ * °^ s Pi^ es °f neuron j in the k ISI of neuron i. 

^k 

(60) 

Let us define the firing rate fj^\ of neuron j'(^ i) in the k th ISI of neuron i, and /|*' = A. We 
obtain 

— H (i) , = - /«. f« . (61) 

rp 33' rj, ^ j k J k,j ''k.J * ' 

where T is the duration of the recording. The right hand side of the above equation can be 
interpreted as the covariance matrix of the rates fj, *. , where each ISI of neuron i is weighted 
proportionally to its duration. For vanishing couplings, these instantaneous rates are decoupled 
from neuron to neuron. Hence, fj?\ fluctuates around the average firing rate fj (number of the 
spikes fired by neuron j, divided by T), with a variance we denote by (fj)c- This statement 

holds for j ^ i; in addition we define fi = —. Neglecting terms of the order of r~ 2 , we end up 
with the following approximation for the Hessian matrix, 

°" 2 rx(0 t t , x u \ if 3 ¥= i . / co n 

— H . ., = fi Ui + o,- ,/ Mi where u)j = < 3' .„ . . , (62) 
T H JJ ' J JJ 1 J \ if ] = i . K ' 

which becomes exact in the limit of infinitely long recordings. 

The matrix HW is the sum of a rank one matrix plus a diagonal matrix. For small values 
of <t, the fluctuations of the firing rates, represented by the uij's, are expected to be small 
compared to the product of any two average firing rates. We immediately deduce that the 
largest eigenvector of the matrix, \ m ax, has components (y m ax)j = fj for all j = 1, . . . , N. 
The associated eigenvalue, X m ax, is given by 

f! Xmax = | £(/,) 2 , (63) 
3 

where S is the total number of spikes. If the neurons have quantitatively similar firing rates 
~ (/), then ff 2 A maI /S ~ (/}. The probability density of vector j/ 1 ' is 



P({v (l) }\T) ~ P({e (l) }|T) x 



cxp 



1 y f^)_ e w\ H w Y„(*)_ fi (^ 

2 ^ \ 3 3 ) ],]' \ 3 3' j 



1,3,3 



(64) 
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A fluctuation <5v = {<5(/;t), SJij} around the most likely values for the current and the cou- 
plings, and along vector v max , will change the log-likelihood by 

Jto&P\ (<5v ■ v max ) 2 T ( ^ , ,V (<^ e ) 2 

5 {— )= 2S(v max) * = -2^s[ SIt+ ^ ) Jtifi ) *-2*m- (65) 

Hence the effective current If is associated to the largest eigenmode, and is the parameter 
requiring the least number of data to be inferred. 

We now look for the smallest eigenvalue, A m ; n . Numerical investigations suggest that 
the associated eigenvector, v m i n , correspond to fluctuations of the current Ii only. We thus 
assume that the components of v m i„ are: (v m in)i — 1> an( i (v m in).j = ~ e j with tj -C 1. The 
eigensystem we need to solve is 

°- 2 A min = — ( — - ^2 e j'fj'l ( 66 ) 

T V r j>&i) ) 

-a 2 \ min e 3 =T fj (-- J2 e r fj')- T ^ V i¥ • ( 67 ) 

According to 16711 and (16611 . we have, for all j), 

£ . = /j r CT 2 \ mm + o((t2 ^ ^ (6g) 

Inserting this expression for the components tj of the eigenvector into (1661) . we obtain 



T — 



T 



A min = t -pr — ■ (69) 

1 + £ — i 5 

If all neurons have quantitatively similar firing rates, (/}, and variances, (/ 2 } c , we obtain 
0" 2 ^min/5 ~ (/ 2 )c/(-/V 2 (/} 3 T 2 ). According to 1681 1. the components ej of the eigenvector are 
very small, tj ~ l/(N 2 r(f)), for all j ^ i. Hence v m ;„ is localized on its current component 
only. A fluctuation <5v = {8(Iir), SJij}, where the SJij's are chosen to be orthogonal to all the 
other eigenmodes of H^' , modifies the log-likelihood by 



A mi „ (5v • v mln ) 2 </ 2 >c 



2 5(v mm )2 2<x2Ar2(/>3 



(70) 



We conclude that the current Ii is the hardest parameter to infer, i.e. the one requiring the 
largest number of data. 

When the membrane leaking time becomes of the order of, or smaller than the average 
ISI duration, the above calculation has to be modified. From a qualitative point of view, the 
average firing rate fj must now be defined as the mean number of spikes emitted by the neuron 
j in a time-window of duration r preceding a spike of neuron i, divided by t, see 1201 1. The 
eigenvector of HW with largest eigenvalue \ m ax is still given by (y ma x)j = fj, with /j = 1/r, 
and 

°j A — - jf E(/j) 2 • ( 71 ) 

i 

Again, these fluctuations are associated to the effective current, with the newly defined average 
firing rates /,. As r gets smaller and smaller, all the rates fj with j ^ i become smaller and 
smaller compared to /;, and the effective current If gets closer and closer to the true current 
Ii. Obviously, the inference of the synaptic coupling Jij is possible if the firing rate fj defined 
on a time- window of duration r preceding a spike of neuron i is much larger than 1 /T, 
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D Fluctuations of the potential around the optimal path at small noise 



In this Appendix, we derive formula J 1 41) for the fluctuations of the potential around its optimal 
value at the mid-point of the ISI. A useful formulation for Pfpt m © can be given in terms 
of a path integral over the potential, 



PFPT(ti,k+l I *i,fci {tj,l}, {Jij}> h) — 

d fVi(t- k + 1 )<V th 



(72) 



at 



-/ 

.fc+i J\ 



+ 1 Jv ^ t Zk*>=° 



VVi(t) exp 



2a 2 



ClVi{t);k,T,J,I] 



The measure T>Vi(t) in the path-integral H72II is restricted to the potentials Vi(t) remaining 
smaller than the threshold Vth at all times t. The upper bound Vi(tJ k ) < V t h means that 
the integral is performed over all the values of the potential smaller than Vth a * time 
while Vi(ti k) is constrained to be zero. 

We introduce the dimensionless variable i/>i (t) = (Vi (t) — V* (t))/V t h to represent the time- 
dependent fluctuation of the potential (Fig. [2p). According to (IT2I) and {I}, the log probability 
density of a path- fluctuation ipi(t) on the inter-spike interval [t; ti k+l] i s > after multiplication 
by<x 2 , 



C[A(t)V th + V*(t);k,T,J,X] 

V?, r*i,h+i r 
C[V*(t);k,T,J,X\ + -^ J ' dt' J 



dtipi(t') 



8 2 C 



5V*{t')SV*(t) 



«i(<) + 0(^) 



L*(T\J,1) 



2 



dtipiit) 



d 2 



(73) 



up to an additive term independent of ipi. Note that we have used the optimality condition 
1B1 to exclude terms linear in i/>; in II73I I. We now want to perform the path integral over the 
fluctuations tpi(t) in H72II . When cr is small we may discard the cubic and higher order terms 
in if)i. The boundary condition on ipi are j/)j(tj j.) = fc+l) = 0: the values of the potential 

Vj(t) are constrained right after and before the emission of a spike, and, hence, cannot fluctuate 
(Fig.[2p). We therefore write the fluctuations ipi(t) as the following Fourier series, 



nir(t — tj t k) 



(74) 



where the tp n are stochastic coefficients. The integral on the last line of d 73 1) can be calculated 
with the result 



dt iPi(t) 



d?_ 
dt 2 



p{CV th f 
4r 



(75) 



where 



(76) 



is the duration of the ISI measured in units of the membrane leaking time. Hence, if we relax 
the constraint that the fluctuating potential should remain below threshold at all times, the 
ipn's are independent Gaussian variables with zero means and variances 



2ra 2 



2a 2 p 



(CV th ) 2 p 2 + n 2 n 2 p 2 + n 2 n 2 



(77) 



where a is defined in 113> . We may now calculate the variance of ipi at the mid-point of the 
ISI, see Fig. EE, 



ti,k + *i,fc+l 



2 ^" sin ("T 

n>l 



2J -^2p+l 

p>0 



(78) 



Summing up the series over p in JTSI gives expression J 141 1 . 
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Fig. 11 Cost-energy function U over the current as a function of the ratio I /(gV t h) f° r different 
values of a (A), g (B), and the inter-spike interval St (C). Values of the parameters are: A. 
St/r = .025, and cj/(V th yfgC) = .016, .16, .32, .64 from right to left; B. g = 1,5,10,40 C/St 
from top to down on the left side, with a/\ r 8i = I; C. v/(V th v r gC) = .32, and St/r = 
1, 10, 50, 100 from bottom to up. 



E Expression of the moving threshold and alternative procedures 

In Section 12.61 we explain that the value of the moving threshold, VJi? , is estimated from the 



intersection of the tangent to the probability of survival in V = Vth with the p B 
Hence. 



v t f = v th + U d ^{&t\v = v th ) 



\ line. 



(79) 



The slope of p s can be expressed in terms of a series of parabolic cylinder functions (Alili, 
Patie and Pcrderscn, 2005), 



^(St\V = V th ) 
av 



E exp(-7ij St/r) , 
Hi U 



a \g 



(80) 



where D' n (z) denotes the derivative of Weber's function of order n, D n (z), with respect to its 
argument z. The normalization coefficients are 



dV D n 



(81) 



The orders nj, i = 0, 1, 2, . . ., are the roots of the equation (Mei and Lee, 1983) 



(82) 



with < no < ni < ri2 < • • ■■ The gap between successive levels, rti+i — ni, is larger than 1. 
Note that the contributions from high orders ri{ decay exponentially with St/r in (1806 - Hence, 
in practice, the summation can be carried out over a finite number of terms. 

The Moving Threshold procedure was designed to take into account the effects of a mod- 
erate noise level, <r. An alternative approximate procedure consists in subtracting to the log- 
likelihood a cost-function preventing the current, or the effective current from getting too close 
to gVth - F° r a quantitative treatment consider a single neuron in the absence of synaptic input, 
for which Pfpt can be calculated under the form of a series of parabolic cylinder functions, 
see above. We denote by Pp PT the approximation to Pfpt obtained when taking into account 
the optimal path only. We define the cost-energy function 



U{I;g,(T,T) = log 



p FPT (5t;g,a, I) 



pf pT (5t;g,a,I) 



(83) 
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for the current /. We show in Fig. 1111 the shape of U for different values of g, a, and the 
inter-spike interval St. As expected from above, this cost function is essentially flat when 
I/(gVth) <S 1, and is repulsive when I/{gVth) ~~ ► 1- The repulsion is strong when the inter- 
spike interval, St, the membrane conductance, g, and the noise standard deviation, a, are 
large. 

In presence of synaptic inputs, we approximate the non-perturbative corrections by sub- 
tracting (TV; — 1)1/ (If) to our log-likelihood, where TV; is the number of spikes of neuron i, and 
7? its effective current. This simple approximation preserves the concavity of the log-likelihood 
and is computationally simple since U has to be calculated only once for each step and neuron. 
Simulations show that the performance of the inference algorithm with the cost function U is 
quantitatively similar to the one obtained with the Moving Threshold procedure. 
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